## 
rm(list = ls())

library(tidyverse)
library(lfe)
library(pbapply)
library(broom)

## Parameters

bw_samp = 10000

## Get the data

census <- read_rds('data/muni_finance.rds')%>% 
  mutate(state_id = factor(state_id)) %>% 
  mutate(east = ifelse(as.numeric(substr(ags, 1, 2)) > 11, 1, 0))

## Rename states
levels(census$state_id) <- c('S-H', 'HH', 'NDS', 'Bremen', 'NRW',
                             'Hessen', 'RP', 'BW', 'Bayern', 'Saarland',
                             'Berlin', 'Brandenburg', 'M-V', 'Sachsen',
                             'Sachsen-Anhalt', 'Thueringen')


## To long (for employment)
df_use <- census %>% 
  dplyr::select(pop_dec_09, ags, matches('nat_ft'),
                treated, state_id) %>% 
  pivot_longer(cols = matches('nat_ft'), 
               names_sep = '_(?!.*_)',
               names_to = c('name', 'year')) %>% 
  pivot_wider(names_from = 'name', values_from = 'value') %>% 
  mutate(year = as.numeric(year)) %>%
  filter(between(pop_dec_09, 0, 20000)) %>% 
  filter(!is.na(emp_nat_ft_pc)) %>%  
  mutate(emp_nat_ft_pc = emp_nat_ft_pc * 1000)

## to long (for wages)

df_use_wages <- census %>% 
  dplyr::select(pop_dec_09, ags, matches('destatis'),
                treated, state_id) %>% 
  pivot_longer(cols = matches('destatis_'), 
               names_sep = '_(?!.*_)',
               names_to = c('name', 'year')) %>% 
  pivot_wider(names_from = 'name', values_from = 'value') %>% 
  mutate(year = as.numeric(year)) %>%
  filter(between(pop_dec_09, 10000 - bw_samp, 10000 + bw_samp)) %>% 
  filter(!is.na(destatis_wage))

## Merge

df_use <- df_use %>% 
  left_join(df_use_wages %>% dplyr::select(ags, year, destatis_wage))

## Post and t

df_use <- df_use %>%
  mutate(post = ifelse(year > 2013, 1, 0),
         treated_post = (as.numeric(treated) - 1) * post)

## Get info on census application and timing

states <- read_rds("data/states_census.rds") %>%
  dplyr::select(state_name, applies_census, census_first_year) %>%
  mutate(state_name = dplyr::recode(state_name, `Schleswig-Holstein` = 'S-H',
                                    `Niedersachsen` = 'NDS',
                                    `Nordrhein-Westfalen` = 'NRW',
                                    `Thüringen` = 'Thueringen',
                                    `Baden-Württemberg` = 'BW',
                                    `Mecklenburg-Vorpommern` = 'M-V',
                                    `Rheinland-Pfalz` = 'RP'))

## Merge to DiD file

df_use <- left_join(df_use, states, 
                    by = c('state_id' = 'state_name')) 

## Time relative to treatment

df_use <- df_use %>%
  filter(!is.na(year)) %>% 
  mutate(time_rel = year - census_first_year,
         time_rel = ifelse(time_rel > -1, time_rel + 1, time_rel),
         post_rel = ifelse(time_rel > -1, 1, 0),
         treated_post_rel = post_rel * (as.numeric(treated) - 1))

## Declare outcomes

outcomes <- c("emp_nat_ft_pc", "destatis_wage")

## Iterate over years pre/post census, as well as bandwidths

periods = 1:3
bw_vec <- seq(500, 2500, 500)

out_sens_p <- pblapply(bw_vec, function(i) {
  lapply(outcomes, function(o) {
    lapply(periods, function(p) {
      
      df_temp <- df_use %>% 
        filter(!(state_id == 'BW'))
      
      ## Temp outcome var
      
      df_temp[, "o"] <- df_temp[, o]
      
      ## Run
      
      m4 <- felm(o ~ treated_post_rel | ags + time_rel | 0 | ags,
                 data = df_temp %>% filter(between(pop_dec_09, 10000 - i,
                                                   10000 + i)) %>%
                   filter(!state_id == 'B-W'), 
                 subset = applies_census == 1 & 
                   between(time_rel, -p,
                           p))
      
      ## Get coef
      out <- broom::tidy(m4, conf.int = T) %>%
        mutate(bw = i, period = p,
               outcome = o)
      out
    }) %>% reduce(rbind)
  }) %>% reduce(rbind)
}) %>% reduce(rbind) %>% 
  mutate(conf.low90 = estimate  - qnorm(0.95) * std.error,
         conf.high90 = estimate  + qnorm(0.95) * std.error)

## Prep for plot

pd <- position_dodge(0.4)

## Plot for paper 

plot_paper <- out_sens_p %>% 
  mutate(period = paste0('Years included\npre/post treatment: ', period)) %>% 
  filter(outcome == "emp_nat_ft_pc")

## Figure A.18: effect on full-time employment ----

p1 <- ggplot(plot_paper, 
             aes(bw, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.3, position = pd) +
  geom_errorbar(aes(ymin = conf.low90,
                    ymax = conf.high90),
                width = 0, size = 0.8,
                position = pd) +
  geom_point(shape = 21, fill = 'white', size = 2,
             position = pd) +
  theme_bw() + xlab('BW around 10,000 inhabitants') +
  ylab('Effect on per 1,000 capita\nfull time employment') +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  facet_wrap(~period, ncol = 3) + 
  scale_x_continuous(breaks = unique(plot_paper$bw))
p1 

## Figure A.20: Effect on wages ----

plot_paper <- out_sens_p %>% 
  mutate(period = paste0('Years included\npre/post treatment: ', period)) %>% 
  filter(outcome != "emp_nat_ft_pc")

p2 <- ggplot(plot_paper, 
             aes(bw, estimate)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0, size = 0.3, position = pd) +
  geom_errorbar(aes(ymin = conf.low90,
                    ymax = conf.high90),
                width = 0, size = 0.8,
                position = pd) +
  geom_point(shape = 21, fill = 'white', size = 2,
             position = pd) +
  theme_bw() + xlab('BW around 10,000 inhabitants') +
  ylab('Effect on average monthly wages') +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  facet_wrap(~period, ncol = 3) + 
  scale_x_continuous(breaks = unique(plot_paper$bw))
p2 
